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Abstract 

In this paper, we focus on a spatial Holling-type IV predator-prey model which contains some 
important factors, such as diffusion, noise (random fluctuations) and external periodic forcing. By 
a brief stability and bifurcation analysis, we arrive at the Hopf and Turing bifurcation surface and 
derive the symbolic conditions for Hopf and Turing bifurcation in the spatial domain. Based on 
the stability and bifurcation analysis, we obtain spiral pattern formation via numerical simulation. 
Additionally, we study the model with colored noise and external periodic forcing. From the 
numerical results, we know that noise or external periodic forcing can induce instability and enhance 
the oscillation of the species, and resonant response. Our results show that modeling by reaction- 
diffusion equations is an appropriate tool for investigating fundamental mechanisms of complex 
spatiotemporal dynamics. 
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I. INTRODUCTION 



Predation, a complex natural phenomena, exists widely in the world, e.g., the sea, the 
plain, the forest, the desert and so on |7j. To model this phenomenon, predator-prey model 



has been suggested for a long time since the pioneer work of Lotka and Volterra 



221. The 



predator-prey model is a kind of "pursuit and evasion" system in which the prey try to 



evade the predator and the predator try to catch the prey if they interact [36]. Pursuit 
means the predator try to shorten the spatial distance between the predator and the prey, 
while evasion means the prey try to widen this spatial distance. In fact, the predator-prey 
model is a mathematical method to approximate some part of our real world. And the 
dynamic behavior of the predator-prey model has long been and will continue to be one of 
the dominant themes in both ecology and mathematical ecology due to its universal existence 
and importance jf], [2^ 



Generally, a classical predator-prey model can be written as the form 



Q: 



dN 
dt 



Nf(N)-mPg{N,P), 



dP 
dt 



P[cmg(N,P) -d), 



(1) 



where N and P stand for prey and predator quantity, respectively, f(N) the per capita 
rate of increase of the prey in absence of predation, d the food-independent death rate of 
predator, g(N, P) the functional response, the prey consumption rate by an average single 
predator, which obviously increases with the prey consumption rate and can be influenced by 
the predator density, and which refers to the change in the density of prey attached per unit 
time per predator as the prey density changes, mg(N, P) the amount of prey consumed per 
predator per unit time, and cmg(N, P) the predator production per capita with predation. 
In population dynamics, a functional response g(N,P) describes the relationship be- 



tween the predator and their prey, and the predator-prey mode 
corresponding functional response for its key position 



tion ecology, both ecologists and mathematicians 
predator-prey models [l, 



type I-III, originally due to Holling |18|, 




is always named after the 
42j. In the history of popula- 
lave a great interest in the Holling-type 



55 



56 



including Holling 
19( , and Holling type IV, suggested by Andrews Q . 
The Holling-type functional responses are so-called "prey- dependent" type [lj, for g(N,P) 
in Eq.([T]) is a function only related to prey N. The classical expression of Holling-type II 



functional response is g(N, P) 



mN 
1+bN 



and g(N, P) 



mN 2 
1+aN 2 



is called Holling-type III. The 
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Holling-type IV functional response is written as follows: 



mN , . 



Fuction (j2J) is called Monod-Haldane-type functional response too 



381 ] . In addition, when 



b = 0, a simplified form g(N,P) = 1 ™^f r , is proposed by Sokol and Howell |45|], and some 



scholars also called it as Holling-type- IV 38|, |42|]. In this paper, we focus on the Holling-type 



IV functional response which takes the form (j2J), and the corresponding predator-prey model 
which takes the form: 

dN- rN (l_N\ mNP dP-p(_ a+ cmN \ (3) 

dt ' iV ^ K) 1+bN+aN 2 ' dt r \ 1 ^ 1+bN+aN 2 ) i V ' 

where r > stands for maximum per capita growth rate of the prey, m > the capture 
rate, c > the conversion rate of prey captured by predator, q > the food-independent 
death rate of predator and K > the carrying capacity, a > the so-called half-saturation 
constant, b > —2^fa such that the denominator of above system does not vanish for non- 
negative N. 

On the other hand, the real world we lived in is a spatial world, and spatial patterns 
are ubiquitous in nature, which modify the temporal dynamics and stability properties of 
population density on range of spatial scales, whose effects must be incorporated in temporal 
ecological models that do not represent space explicitly [37]. And the spatial component of 
ecological interactions has been identified as an important factor in how ecological commu- 
nities are shaped. Empirical evidence suggests that the spatial scale and structure of the 
environment can influence population interactions and the composition of communities [3] • 

Reaction-diffusion model is a typical spatial extended model. It involves not only time 
but also space and consists of several species which react with each other and diffuse within 
the spatial domain. It also involves a pair of partial differential equations, and represents the 
time course of reacting and diffusing process. In the spatial extended predator-prey model, 
the interaction between the predator and the prey is the reaction item, and the diffusion 
item comes into being for the predator's "pursuit" and the prey's "evasion" . Diffusion is a 
spatial process, and the whole model describes the evolution of the predator and the prey 
going with time. 



Decades after Turing 48] demonstrated that spatial patterns could arise from the inter- 



action of reactions or growth processes and diffusion, react ion- diffusion models have been 



studied in ecology to describe the population dynamics of predator-prey model for a long 
time since Segel and Jackson applied Turing's idea 43]. Since then, a new field of ecology, 
pattern formation, came into being. The problem of pattern and scale is the central prob- 
lem in ecology, unifying population biology and ecosystem science, and marrying basic and 
applied ecology [2^] . And the study of spatial patterns in the distribution of organi sms is a 



central issue in ec olog 
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geo 
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chemistry. 
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physics and so on 
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13, 



571 ] . Theoretical work has 



shown that spatial and temporal pattern formation can play a very important role in ecolog- 
ical and evolutionary systems. Patterns can affect, for example, stability of ecosystems, the 
coexistence of species, invasion of mutants and chaos. Moreover, the patterns themselves 
may interact, leading to selection on the level of patterns, interlocking ecoevolutionary time 
scales, evolutionary stagnation and diversity. 

Based on the above discussions, the spatial extended Holling-type IV predator-prey model 
with reaction diffusion takes the form: 



8N — r Kf (1 _ K) mm , a, v 2 N 

dt - riV \ L K) 1+bN+aN? +«1 V iV ' 



dP 
dt 



P(-q 



cmN 



l+bN+ 



^) + d 2 V 2 P, (4) 



where d\ and 62 are the diffusion coefficients respectively, V 2 = ^ + ^2 the usual Laplacian 
operator in two-dimensional space, and other parameters are the same definition as those in 
model p. 

Easy to know that, when a spatial extended predator-prey model is considered, the evo- 
lution of the model is decided by two sorts of sources (internal source and external source) 
which act together. The internal source is the dynamics of the individuals of the model, and 
the external source is the variability of environment. Some of the variability is periodic, such 
as temperature, water, food supply of the prey and mating habits. It is necessary and impor- 
tant to consider models with periodic ecological parameters or perturbations which might 



be quite naturally exposed 



101 ] . These periodic factors are regarded as the external peri- 



odic forcing in the predator-prey systems. The external forcing can affect the population of 
the predator and prey, respectively, which would go extinct in a deterministic environment. 
And some of the variability is irregular, such as the the seasonal changes of the weather, 
food supply of the prey, mating habits, and the affects of this variability are called "noise" . 
Ecological population dynamics are inevitably "noisy" 22j. In the predator-prey systems, 
the random fluctuations are also undeniably arising from either environmental variability or 
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internal species. To quantify the relationship between fluctuations and species' concentra- 
tion with spatial degrees of freedom, the consideration of these fluctuations supposes to deal 
with noisy quantities whose variance might at time be a sizable fraction of their mean levels. 
For example, the birth and death processes of individuals are intrinsically stochastic fluctu- 



ations which becomes especially pronounced when the number of individuals is small [34]. 
Moreover, there are many other stochastically factors causing predator-prey population to 
change, such as effects of spatial structure of the habitat on the predator-prey ecosystem. 
The interaction between the predator and prey, which are far from being uniformly dis- 
tributed, also introduce randomness. And these processes can be regarded as parameter 
fluctuates irregularly with spaces and time. 

The induced effects of the external forcing and noise in population dynamics, such as 
pattern formation, stochastic resonance, delayed extinction, enhanced stability, quasi peri- 



odic oscillations and so on, have 
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3een investigated with an increasing interest in the past 
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571 ] . And noise cannot systematically be neglected 



decades _ 

in models of population dynamics [50]. Zhou and Kurths [57] concluded these periodic vari- 
abilities as external forcing, and investigated the interplay among noise, excitability, mixing 
and external forcing in excitable media advected by a chaotic flow, in a two-dimensional 
FitzHugh-Nagumo model described by a set of reaction-advection-diffusion equations. And 
Si et al [44 [studied the propagation of traveling waves in sub-excitable systems driven, and 
Liu et al [28( considered a spatial extended phytoplankton-zooplankton system with addi- 
tive noise and periodic forcing. Following these models they considered, the Holling-type IV 
predator-prey model with external periodic forcing and colored noise is as follows: 

f = rN (1 - £ ) - + A sin (ut) + d^N, 

% = P(~1+ ITWW) + V(r, t) + d 2 V 2 P, 

where A sin(a>£) denotes the periodic forcing with amplitude A and angular frequency u. The 
colored noise term rj(r, t) ( r = (x, y) ) is introduced additively in space and time, referring to 
the fluctuations in the predator increase rate, which partially results from the environmental 
factors such as epidemics, weather and nature disasters, and it is the Ornstein-Uhlenbech 
process that obeys the following linear stochastic partial differential equation 

«2M = -i,( r , t ) + i {(r , t ), (6, 
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where £(r, t) is a Gaussian white noise or the so called Markovian random telegraph process 
in both space and time with zero mean and correlation: 

(£(r, t )) = 0, (£(r, t)£(r', t')) = 2e5(r - v')S(t - t'), 

where (•) denotes mean value with respect to the noise £(r, t), and 5 the Dirac delta-function, 

5(r — r') the spatial correlation function of the Gaussian white noise £(r, t) . 

Integrating equation (jBJ) with respect to time t, we get 

t 1 t f l 
rj(r, t) = r)(r, 0)e~~ H — / e~^(r,s)ds. 

T Jo 

The mean value of the colored noise is 

1 /"* 

( V (r,t)) = (r,(r,0))e" + -e~$ / e^(£(r, s))ds = ( V (t, 0))e^, 

T Jo 

and the correlation function of the colored noise is given by 

( V (v,t) V (r',t')) = ( 7? ( r ,0))(r ? (r',0))e-^ + ^e-^ : £^ e^(£(r, S )£(r', s'))dsds' (7) 

= (r/(r,0))(r/(r / ,0))e-^ + ^e-^5(r-r') f f e^8(t - t')dsds r (8) 

T Jo Jo 

= (r){r, 0)) (77(r', 0))e"^ + -8{v - r / )(e" 1 ^ - 2e~^ + e"^). (9) 

r 

Let t — > +oo, then 

(?7(r, t)r](r', t)) — > — e~^^5(r — r ). 
r 

The colored noise f](v,t) generated in this way represents a simple spatiotemporal struc- 
tured noise that can be used in real mimic situations, which is temporally correlated and 
white in space, and satisfies 

(r,(v,t)r](r',t')) = -e-^6(r - v'), 

T 

where the temporal memory of the stochastic process controlled by r and e is the intensity 
of noise. In this paper, we set r = 1. 

Based on these discussion above, in this paper, we mainly focus on the spatiotemporal 
dynamics of model Q and (J5j). And the organization is as follows: In section 2, we employ 
the method of stability analysis to derive the symbolic conditions for Hopf and Turing 
bifurcation in the spatial domain. In section 3, we give the complex dynamics of model (jlj) 
and (jnj), involving pattern formation, phase portraits, time series plots and resonant response 
and so on, via numerical simulation. Then, in the last section, we give some discussions and 
remarks. 
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II. HOPF AND TURING BIFURCATION 



The non-spatial model ([3]) has at least two equilibria (steady states) which correspond to 
spatially homogeneous equilibria of the model (jl]) and ([5]), in the positive quadrant: (0,0) 
(total extinct) is a saddle; (K, 0) (extinct of the predator, or prey-only) is a attracting node 
if Q > i+^b+aK^ » a saddle if 1 < i+xl+aK^ or a saddle-node if q = 1+ ^ aR2 - When 

(a, b, c, m, q, r, K) £ E\, here, £i = {(a, 6, c, m, q, r, K) \ mc > qb, q 2 a < (mc — qb) 2 , 

^(mc-qb) 2 -Aq 2 a > ~ m2<a + 2 qbmc tZ?+™b+faf b+2 > ~ ^ + K < °h 
there exists unique stationary coexistent state (Ni*,Pi*), where 

at* 1 — qb+mc— A p * cr((—mc+bq+qaK )Ni*+q) 

iVl ~~ 2 ~q~a ' r l ~~ ^pi? • 

On the other hand, when 
(a, b, c, m, q, r, if) £ P 2 , here, E 2 = {(a, 6, c, m, g, r, K)| mc > qb, q 2 a < (mc — qb) 2 , 

J(mC-qb) 2 -Aq 2 a > - ~ m ^^ 2 gbrnc+qaKrnc-q^aKb+2 gWjg > 0, k _ vac + R q| 
V V ^ / 3 —mc+qb+qaK 'a ga J ' 

there exists another unique stationary coexistent state (N 2 *,P 2 *) implying: 

at * 1 —qb+mc+A r> * er((— me +bq+qaK )A r 2*+g) 

iV 2 - 2 ^ ' ^2 - ^ • 

It is worth mentioning that equilibria (N^,P*) and (N£,P£) cannot coexist. In this 
paper, we mainly focus on the dynamics of (N£, P*) and rewrite it as (N*, P*). The dynamic 
behavior of (N 2 , P 2 * ) is similar to those of (N*, P±). 

To perform a linear stability analysis, we linearize model (j3]) around the stationary state 
(N* , P*) for small space- and time- dependent fluctuations and expand them in Fourier space 

N(r, t) ~ N*e xt e llr , P(v, t) ~ P V f e* r , r = (ar, y), k = (k x , k y ). 

where A is the eigenvalue of the Jacobian matrix of model ([3]). 

Hopf bifurcation is an instability induced by the transformation of the stability of a focus. 
Mathematically speaking, Hopf bifurcation occurs when Im(A) ^ 0, Re(A) = at k — 0, 
where Im(A) is the imaginary part, Re(A) the real part and k the wave number. So we get 
the Hopf bifurcation surface 

H = {(a,b,c,m,q,r,K)\ det(J ) > 0,trace(J ) = 0}, (10) 
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where 

Hp+fn- (t 9rN*\ n | mqP*+cm(r-2^)N* mq N* P* (b+2 aN*) 

det(J ) - - [r - 2 — ) q + —^———^ ; 

. ^ rAr , m(-P*+cAT*+aV* 2 P*+bc7V* 2 +eiV* 3 a) 

trace( J ) = r — 2 ^ g H - — 



(l+bN*+aN* 2 ) 

the frequency of periodic oscillations in time uo H satisfies uo H = Im(A) = y^det(J Q ), and the 
corresponding wavelength A# satisfies \h = = ^//^j ) ' Especially, we ^ a ke K as the 



bifurcation parameter, and can get the critical value of Hopf bifurcation from Eq. f jlOT) : 



Kh = (—(aq 2 (5 mc — 3 qb) — (3 mc — qb) (mc — qb) 2 ) \J (mc — qb) 2 — 4 q 2 a 

—4 q A a 2 + q 2 (mc — qb) (11 mc — 5 gfe)a — (3 mc — qb)(mc — qb) 3 )/ (—aq(((2 mc — g&) 

(mc — qb) — 2 q 2 a) a/ (mc — g&) 2 — 4 g 2 a — 2aq 2 (3mc — 2 qb) + (2mc — qb)(mc — qb) 2 )). 

(11) 

Turing instability is induced only by "pursuit and evasion" if the predator can catch the 
prey by pursuit. We call the critical state of Turing instability as Turing bifurcation. Turing 
bifurcation occurs when Im(A) = 0, Re(A) = at k = fcy 7^ 0, and the wavenumber kr 
satisfies k\ = \J dc ^f ■ In addition, at the Turing threshold, the spatial symmetry of the 
system is broken and the patterns are stationary in time and oscillatory in space with the 
wavelength At = And the Turing bifurcation surface implies: 

T= {(a, b, c, m, q, r, di,cf 2 ,-^)| det( J k ) = 0, trace( J k ) = 0}, (12) 

where 

det(J fe ) = - (r - 2 rfl - d x k 2 ) (q + d 2 k 2 ) ^^Ktbtlf 

m(b+2 aN*)(q+d 2 k 2 )N* P* 
{l+bN*+aN* 2 ) 2 ' 

m(-P*+cN*+aN* 2 P*+bcN* 2 +cN* 3 a) 



trace( J k ) = r - 2 ^ - q - (d x + d 2 ) k 2 + 



K ' -'" (l+bN*+aN* 2 Y 



the critical value of Turing bifurcation can be obtained from Eq. (fT2]) as follows 



where 



Kt = (13) 



Fi = r((4q 2 a(2mc — qb) — (3mc — qb)(mc — qb) 2 ) \J (mc — qb) 2 — 4 q 2 a + ((mc — qb) 2 — 4q 2 a) 
• ( (3 mc — qb) (mc — qb) — 2 q 2 a) ) ( (3 mc — qb) (mc — qb) — Aq 2 a 
— (3mc + g&) a/ (mc — g&) 2 — 4 g 2 a )<i 2 , 
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F 2 = qa(2 mc((mc — qb) \J (mc — qb) 2 — 4 q 2 a + 4 q 2 a — (mc — qb) 2 )B + ((3 mc — qb)d 2 
■(2 m 2 c 2 — 3 qbmc + q 2 b 2 — 4 g 2 a)r + 2 qd\mc(mc — qb) 2 )((mc — qb) 2 — 4 g 2 a) 
+(-2d 2 ((3mc — qb)(2mc — qb)(mc — qb) 2 — 2q 2 a(—Aq 2 a + 3q 2 b 2 — 12 qbmc 
+ 11 m 2 c 2 ))r — 4qcm(mc — qb)d\((mc — gfe) 2 — 4 q 2 a))\J (mc — g&) 2 — 4 g 2 a 
+ (g 2 6 2 — 2 gfomc + m 2 c 2 — 4q 2 a)((2mc — qb)d 2 (3m 2 c 2 — A qbmc + q 2 b 2 — 4q 2 a)r 
+2mcqdi((mc — qb) 2 — 4g 2 a))), 

B = (— 2 diq(((d\bq 2 — qd\mc — qrbd 2 + 2rmcd2)( / lq 2 a — (mc — qb) 2 ) — rmcd^imc — qb) 2 ) 

■ (mc — qb) 2 — 4 q 2 a + (8 q^d 2 r — 8 g 5 <i 1 )a 2 + 4 (mc — qb)aq 2 rmcd2 + (3 rmcd 2 

—qd\mc + c?i6g 2 — qrbd 2 )((mc — g6) 3 — 6 (mc — qb)aq 2 )))^ . 
Linear stability analysis yields the bifurcation diagram with r = 1, a = 0.125, 6=1, 
c = 0.7, m = 0.625, q = 0.18 and rf 2 = 0.2 shown in Figure [T]( A) . In this case, parame- 
ters (a,b,c,m,q,r, K) G E\, and (N*,P*) is the unique stationary coexistent state. From 
Figure [H(a), one can see that the Hopf bifurcation line and the Turing bifurcation curve 
separate the parametric space into three distinct domains. In domain I, located below all 
two bifurcation lines, the steady state is the only stable solution of the model. Domain II is 
the region of pure Hopf instability. When the parameters correspond to domain III, located 
above all two bifurcation lines, both Hopf and Turing instability occur. Figure [H(b) displays 
the dispersion relations showing unstable Hopf mode, transition of Turing mode from stable 
to unstable for model (T4]), e.g., as di decreased. Figure [T](c) illustrates the relation between 
the real and the imaginary parts of the eigenvalue A with K = 2.8 > Kh = 2.279, located 
in domain II, one can see that when k = 0, Re(\(k)) > and lm(\(k)) ^ 0. Figure []Jd) 
displays the case of the critical value of Turing bifurcation K = Kt = 3.499, in this case, 
Re(A(fc)) = and Im(A(A;)) = at k = k T = 2.080. When K = 4.0, located in domain III, 
figure C](e) indicates that at k = 0, Re(X(k)) > 0, lm(\(k)) ^ 0. 

III. SPATIOTEMPORAL DYNAMICS OF THE MODELS 

In this section, we perform extensive numerical simulations of the spatially extended 
model (jlj) and (Q in two-dimensional spaces, and the qualitative results are shown here. All 
our numerical simulations employ the Zero-flux Neumann boundary conditions with a system 
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(c) (d) (e) 



FIG. 1: (a) K — d\ Bifurcation diagram for model dH) with r = 1, a = 0.125, b = 1, c = 0.7, 
m = 0.625, q = 0.18 and d 2 = 0.2. Hopf and Turing bifurcation lines separate the parameter 
space into three domains. Pi (0.02, 2.8), P 2 (0. 02, 3.499), P 3 (0.02, 4.0) are the different selections 
corresponding to (c), (d), (e), respectively, (b) Dispersion relations showing unstable Hopf mode, 
transition of Turing mode from stable to unstable for the system model (jJJ) as d\ decreased. The 
other parameters in (c)-(e) are: d\ = 0.02, the bifurcation parameter K equals: (c) 2.8 > Kh = 
2.279; (d) 3.499 = K T ) (e) 4.0 > K T > K H . The real parts Re(A) and the imaginary parts Im(A) 
are shown by solid curves and dashed curves, respectively. 

size of 200x200 space units. The parameters are r = 1, a = 0.125, b = 1, c = 0.7, m = 0.625, 
q = 0.18, d x = 0.02, d 2 = 0.2 and K = 2.8 or K = 4.0, which satisfy (a, 6, c, m, q, r, K) G E. 
Model (j3j) and (jSJ) are integrated initially in two-dimensional space from the homogeneous 
steady state, i.e., we start with the unstable uniform solution (N*,P*) with small random 
perturbation superimposed, in each, the initial condition is always a small amplitude random 
perturbation (±5 x 10~ 4 ), using a finite difference approximation for model (j3J) or Fourier 
transform method for model (jSJ) for the spatial derivatives and an explicit Euler method for 
the time integration with a time stepsize of At = 1/24 and space stepsize (lattice constant) 
Ax = Ay = 1. We have taken some snapshots with white (black) corresponding to the high 
(low) value of prey N. 
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In the numerical simulations, different types of dynamics are observed and we have found 
that the distributions of the predator and prey are always of the same type. Consequently, 
we can restrict our analysis of pattern formation to one distribution. In this section, we 
show the distribution of prey N, for instance. 

A. Pattern formation of model (4) 

Figure [2] shows the evolution of the spatial patterns of prey N at t = 0, 100, 300, 500, 
1000, 2000, with random small perturbation of the equilibrium (N*,P*) = (0.748,2.132) of 
model (jl]) with K = 2.8, located in domain II, more than the Hopf bifurcation threshold 
Kh = 2.279 and less than the Turing bifurcation threshold Kt = 3.499. In this case, pure 
Hopf instability occurs. One can see that for model (J4j), the random initial distribution (c.f., 
figure E](a)) leads to the formation of macroscopic spiral patterns (c.f., figures EJ^d) to (f)). 
In other words, in this situation, spatially uniform steady-state predator-prey coexistence 
is no longer. Small random fluctuations will be strongly amplified by diffusion, leading to 
nonuniform population distributions. From the analysis in section 2, we find with these 
parameters in domain II, the spiral pattern arises from Hopf instability. The lower panel 
in figure [2] shows the corresponding (g) time series and (h) phase portraits. Figure EJ^g) 
illustrates the evolution process of prey N, periodic oscillating in time finally, (h) exhibits 
that a limit cycle arises, which is caused by the Hopf bifurcation. 

When K = 4.0 > Kt > Kh, in this case, parameters in domain III (figure [T^A)), 
both Hopf and Turing instabilities occur. The nontrivial stationary state is (N*,P*) = 
(0.748,2.365). As an example, the formation of a regular macroscopic two-dimensional 
spatial pattern is shown in figure [3j The lower panel in figure [3] shows the corresponding (g) 
time series plots and (h) phase portraits. 

Comparing this situation (figure [3]) with the one above (figure [2]), easy to see that the 
pattern formations are all spiral wave. From the analysis in section 2, we know that when 
K = 2.8, the wavelength A = 3.100 while at K = 4.0, A = 3.021. And the frequency of 
periodic oscillations in time is as inverse proportion to wavelength, so we can know that 
Turing instability has a positive effect on the frequency and negative effect on wavelength. 
This is the reason why the spiral curves are more dense in figure [3]T) than in figure [2](f). 
On the other hand, one can see that when K = 4.0, the time series plots (c.f., figure [3(g)) 
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FIG. 2: Grey-scaled snapshots of spatiotemporal pattern of the prey N of model with K = 2.8. 
(a) t = 0, (b) t = 100, (c) t = 300, (d) t = 500, (e) t = 1000, (f) t = 2000. The lower panels show 
the corresponding (g) time series plots and (h) phase portraits. 



indicate that when Turing instability occurs, the solution of model (jlj) is strongly oscillatory 
in time while with K = 2.8 (pure Hopf bifurcation emerges) it is periodic (c.f., figure 121(g))- 
In addition, comparing figure [21(g) with figure M^g) , one can see that Turing instability has a 
positive effects on the amplitude of prey N. And from figure [3](h) , one can see that a quasi 
limit cycle emerges while in figure[2](h), it is a cycle. Although there are some different points 
between figure [2] and figure [31 but we can know that Turing instability can't give birth to 
different type patterns. In our previous work 5l|] , we find that Turing instability can change 
pattern type. This may be a very point between the Holling-type IV with Michael-Menton 
functional response of predator-prey model. 

On the other hand, the basic idea of diffusion- driven instability in a reaction-diffusion sys- 
tem can be understood in terms of an activator- inhibitor system or predator- prey model (jlj). 
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FIG. 3: Grey-scaled snapshots of spatiotemporal pattern of the prey N of model with K = 4.0. 
(a) t = 0, (b) t = 100, (c) t = 300, (d) t = 500, (e) t = 1000, (f) t = 2000. The lower panels show 
the corresponding (g) time series plots and (h) phase portraits. 

The functioning of this mechanism is based on three points First, a random increase of 
activator species (prey, N) should have a positive effect on the creation rate of both acti- 
vator (prey, N) and inhibitor (prey, P) species. Second, an increment in inhibitor species 
should have a negative effect on formation rate of both species. Finally, inhibitor species P 
must diffuse faster than activator species N. Certainly, the reaction-diffusion predator-prey 
model PI, with Holling-type IV functional response and predators diffusing faster than prey 
(i.e., d,2 > di), provides this mechanism. And spirals and curves are the most fascinating 
clusters to emerge from the predator-prey model. A spiral will form from a wave front when 
the prey line (which is leading the front) overlaps the pursuing line of predator 17|. The 
prey on the extreme end of the line stops moving as there is no predator in the immediate 
vicinity. However, the prey iV and the predator P in the center of the line continue moving 
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forward. This forms a small trail of prey at one (or both) ends of the front. These prey 
start breeding and the trailing line of prey thickens and attracts the attention of predator 
at the end of the fox line that turn towards this new source of prey. Thus a spiral forms 
with predator P on the inside and prey N on the outside. If the original overlap of prey 
occurs at both ends of the line a double spiral will form. Spirals can also form as a prey 
blob collapses after the predator eats into it. This is the reason why the pattern formation 
of model fllj is spiral wave. 



B. The effect of noise only of model (5) 

Now, we turn our focus on the effect of noise on the predator P of stochastic model (0). 
In this case, A = 0, i.e., the periodic forcing is not at presence. 

Figure 0] shows the dynamics of model ([5]) with noise on the predator. The first row of 
figure HI i.e., (a), e = 0.0001; the second row, (b), e = 0.01; the third row, (c), e = 0.05; and 
the last row of figure HI (d), e = 0.1. And the first collum of figure HI marked as (i), shows 
the snapshots of spatiotemporal pattern of model (jSJ) at t = 2000 with different intensity of 
noise, respectively. In this case, one can see that the pattern formation turns into spatial 
chaotic from spiral wave with the increase of noise intensity e. And the second collum of 
figure HI marked as (ii), displays the phase portraits of model (J3j) with different intensity of 
noise, respectively. We can see that, as noise intensity e increasing, the symmetry of the 
limit cycle is broken and gives rise to chaos. The last collum of figure HJ (hi), illustrates the 
time-series plots of prey iV with different intensity of noise, respectively. One can see that 
noise breaks the periodic oscillations in time and gives rise to drastically ruleless oscillations 
in time. 



C. The effect of periodic forcing of model (5) 



In the previous subsection, we have shown the effect of noise on the predator P of 
model (JSJ). An interesting question is whether such noise-sustained oscillations can be en- 
trained by a weak external forcing, in this case, e — 0, which is investigated here. 



When model ([5 
33, 44. 



response [28 



is noise free, there is a phenomenon of frequency locking or resonant 



571 ] . That is, without noise, the spatially homogeneous oscillation does 
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FIG. 4: Dynamics of model (J!}, for the following noise intensity, (a) e = 0.0001; (b)e = 0.01; 
(c)e = 0.05; (d)e = 0.1. (i) snapshots of pattern formation at time 2000; (ii) phase portraits; (iii) 
time-series plots. ^4 = and the other parameters are the same as those in figure [2j 



not respond to the external periodic forcing when the amplitude A is below a threshold 
whose value depends on the external periodic T in = 2£. Above the threshold, model (jSJ) 
may produce oscillations about period T out with respect to external period T in , which is 
called frequency locking or resonant response. That is, the model produces one spike within 
each of the M = ^fr 1 periods of the external force, called M : 1 resonant response 44], [571 . 

in . ... 



The phenomenon of coherent resonance is of great importance 33j. Following Si 44J, in the 
present paper, the output periodic T out is defined as follows: Tj is the time interval between 
the ith spike and (i + l)th spike, m spikes are taken into account and the average value of 
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them is T out = J2 %/ (m - 1). 

As an example, with the amplitude A = 0.001, figure [5] shows 5 : 1 resonant response with 
uj = 0.27T (a) and u = 0.02n (c), respectively. And figure Mjo) and (d) are the phase portraits 
corresponding to (a) and (c). We can see that when uj = 0.2n, there exists a periodic orbit, 
while uj = 0.027T, a periodic-2 orbit of model © emerges. Obviously, different uo can form 
the same resonant response, and different phase orbits, i.e., different numerical solution of 
model (jSJ), may correspond to the same resonant response. 
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FIG. 5: External periodic forcing induced frequency locking of model(TjI]). The solid curve is time 
series of prey U, the dash curve is the corresponding external periodic forcing. Other parameters 
are the same as those in figure El 



D. The effect of noise and periodic forcing of model (5) 

Now, we consider the dynamics about resonant response of model (jSJ) with both noise 
and periodic forcing. As depicted in figure El the prey can generate 5 : 1 (a), 4 : 1 (c) locked 
oscillations, depending on the amplitude A and angular frequency u. Figures M (b) and 
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(d) illustrate the spiral pattern at t — 2000 corresponding to (a) and (c), respectively. For 
contrast, we change one of the parameters of figure E](c) A = 0.001 to A = 0.01 (e), one can 
see that the resonant response vanishes, the corresponding spiral pattern (f) is similar to 
(b). It indicates that the amplitude A is a control factor for pattern formation. In addition, 
comparing figures E(b) with (d), one can see that the pattern formations are determined by 
noise intensity e, too. 
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FIG. 6: Dynamics of model ([5]) with both noise and periodic forcing. (b)(d)(f) are snapshots at 
t = 2000 corresponding to the left hand side resonant response. The other parameters are the same 
as those in figure [3l 



In Figure [?l we have shown a typical pattern formation process in the 5 : 1 frequency 
locking regime with A = 0.001 and u = 0.2tc. From t = 1870 (a) to t = 1920 (f), the pattern 
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FIG. 7: Typical pattern formation of the forced noisy prey in the 5 : 1 locking region at A = 0.001 
and e = 0.0001 corresponding to figure [6^a). The lower panel shows the time series of the prey iV 
(the solid curve) and the corresponding external periodic forcing (the dash curve) corresponding 
to the snapshots of the patterns. The grey scale, from black to white, is in [0, 3.5] for all the 
snapshots. 

formation of prey N is spiral wave and some small excitations already develop. One can see 
that, during the second periodic of the forcing, the prey is almost fully synchronized and 
relaxes slowly back to the state at moment (f). Obviously, the external periodic forcing at 
moment (e) repeats that at moment (a). However, the prey N does not exactly repeat due 
to a small fluctuation of the phase difference. 

IV. CONCLUSIONS AND REMARKS 

In this paper, we present spatial Holling-type IV predator-prey model containing some 
important factors, such as noise (random fluctuations), the external periodic forcing and 
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diffusion processes. And the numerical simulations were consistent with the predictions 
drawn from the bifurcation analysis, that is, Hopf bifurcation and Turing bifurcation. 

If the parameter K, the carrying capacity, located in the domain II of figure (Ha), the Hopf 
instability occurs, the destruction of the pattern begins from the prey N, while it begins 
from the predator P if K located in the domain III, both Hopf and Turing instabilities 
occur. 

Furthermore, we demonstrate that noise and the external periodic forcing play a key 
role in the predator-prey model (jSJ) with the numerical simulations. We provoke qualitative 
transformations of the response of the model by changing noise intensity, and noise can 
enhance the oscillation of the species density, and format large clusters in the space. The 
periodic oscillations appear when the spatial noise and external periodic forcing are turned 
on. It has also been realized that model (JSJ) is very sensitive to external periodic forcing 
through the natural annual variation of prey growth. In conclusion, we have shown that 
the cooperation between noise and external periodic forcing inherent to the deterministic 
dynamics of periodically driven models gives rise to the appearance of resonant response. 

Significantly, model © exhibits oscillations when both noise and external forcing are 
present. This means that the predator population may be partly due to the external forcing 
and stochastic factors instead of deterministic factors. Therefore, the model for spatially 
extended systems composed by two species could be useful to explain spatiotemporal be- 
havior of populations whose dynamics is strongly affected by noise and the environmental 
physical variables, and the results of this paper are an important step toward providing the 
theoretical biology community with simple practical numerical methods, for investigating 
the key dynamics of realistic predator-prey models. 
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